About the data

Indonesia earthquakes datasets

These datasets contain spatio-temporal information about earthquakes in Indonesia.

quakes_df.rds contains records for individual earthquakes. provinces_df.rds contains records for each province, including aggregations of earthquake counts and densities.

Steps to take

  1. Read in the two data sets
  2. Subset the data to contain only quakes in November 2023, small enough to experiment with
  3. Make a map with earthquake locations plotted
  4. Make an interactive plot of magnitude and date
  5. Link the two plots so that a user can select points in one plot and examine where they are in the other
  6. Develop an animation of earthquakes over time
  7. Create a choropleth map of earthquakes by province
  8. Convert the map to a cartogram
  9. Add interactive mouse-over
  10. Create a quarto document to start a report on earthquakes in Indonesia

1. Read

quakes_sf <- readRDS("data/earthquakes/quakes_df.rds")
provinces_sf <- readRDS("data/earthquakes/provinces_df.rds")

# Check map
ggplot(provinces_sf) + geom_sf()

2. Subset the data

quakes_df <- st_drop_geometry(quakes_sf)
quakes_df <- quakes_df |>
  bind_cols(st_coordinates(quakes_sf)) 
quakes_11_2023 <- quakes_df |>
  filter(year(ymd) == 2023, month(ymd) == 11)

3. Make a map

ggplot(provinces_sf) + 
  geom_sf() +
  geom_point(data=quakes_11_2023, aes(x=X, y=Y)) +
  theme_map()

4. Make an interactive plot

p1 <- ggplot(quakes_11_2023, aes(x=ymd, y=mag)) +
  geom_point()
ggplotly(p1)

Note: sf objects don’t work with plotly if there are multiple geometry types in them.

quakes_11_2023_shared <- highlight_key(quakes_11_2023)

p1 <- ggplot(quakes_11_2023_shared, aes(x=ymd, y=mag)) +
  geom_point() +
  theme_minimal()

gp1 <- ggplotly(p1, width=400, height=400) |>
  config(displayModeBar = FALSE) |>
  highlight(on = "plotly_selected", 
            off = "plotly_doubleclick") 

# multiple types of geometry in provinces_sf needs fixing
provinces_sf2 <- st_cast(provinces_sf, "MULTIPOLYGON")
p2 <- ggplot() + 
  geom_sf(data=provinces_sf2) +
  geom_point(data=quakes_11_2023_shared, aes(x=X, y=Y)) +
  theme_map()

gp2 <- ggplotly(p2, width=800, height=400) |>
  config(displayModeBar = FALSE) |>
  highlight(on = "plotly_selected", 
            off = "plotly_doubleclick") 

bscols(gp1, gp2, widths = c(4, 8))

7. Create a choropleth map

quakes_2023 <- quakes_df |>
  filter(year(ymd) == 2023) |>
  group_by(province) |>
  summarise(mag = mean(mag, na.rm=TRUE),
            n = n(), na.rm=TRUE)

provinces_sf2 <- left_join(provinces_sf2, quakes_2023)

p3 <- ggplot(provinces_sf2) + 
  geom_sf(aes(fill=mag, label=province)) +
  theme_map() +
  scale_fill_viridis_c()

8. Convert to a cartogram

Population data can be found here.

pop <- read_csv("data/province_pop.csv") |>
  rename(province = Provinsi)
provinces_sf2 <- left_join(provinces_sf2, pop)

indonesia_carto <- provinces_sf2 |> 
  st_transform(3395) |> 
  cartogram_cont("2023") |>
  st_transform("WGS84") 

indonesia_carto <- st_cast(indonesia_carto, "MULTIPOLYGON")

ggplot(indonesia_carto) + geom_sf()

indonesia_carto <- left_join(indonesia_carto, quakes_2023)

p4 <- ggplot(indonesia_carto) + 
  geom_sf(aes(fill=mag, label=province)) +
  theme_map() +
  scale_fill_viridis_c()

9. Interactivity

ggplotly(p3)

ggplotly(p4)

10. Embed your best examples into a quarto document

Once you have your quarto document rendering, change the data to use 2022, and re-render, to see how fast it is to re-do the analysis if the data changes.

Resources

Creative Commons License This work is licensed under a Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International License.